m 
O 
o 



O 



Equilibrium Models of Galaxy Clusters with Cooling, Heating 

and Conduction 

M.Briiggen 1 ' 2 

m . brueggenOiu-bremen . de 

ABSTRACT 



It is generally argued that most clusters of galaxies host cooling flows in which 
radiative cooling in the centre causes a slow inflow. However, recent observations 
by Chandra and XMM conflict with the predicted cooling flow rates. Amongst 
other mechanisms, heating by a central active galactic nucleus and thermal con- 
duction have been invoked in order to account for the small mass deposition rates. 
Here, we present a family of hydrostatic models for the intra-cluster medium 
where radiative losses are exactly balanced by thermal conduction and heating 
by a central source. We describe the features of this simple model and fit its 
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parameters to the density and temperature profiles of Hydra A. 



Subject headings: galaxies: active - galaxies: clusters: cooling flows - X-rays: 
galaxies 
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1. Introduction 

The X-ray surface brightness of many clusters of galaxies shows a strong central peak 
which has been interpreted as the signature of a cooling flow (Cowie & Binney 1977, Fabian 
& Nulsen 1977, Sarazin 1988, Fabian 1994). However, the simple cooling flow model conflicts 
with a growing number of observations that show that while the temperature is declining 
in the central region, gas with a temperature below ~1 keV is significantly less abundant 
than predicted. The absence of a cool phase in cooling flows has proven a persistent puzzle. 
Recently, two main candidates for the heating of the gas in the central regions of clusters 
have emerged: (i) heating by outflows from active galactic nuclei (AGN) (Tabor & Binney 
1993, Churazov et al. 2001a, Binney 2001, Briiggen et al. 2002, Briiggen 2003), and (ii) 
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transport of heat from the outer regions of the cluster by thermal conduction (Zakamska & 
Narayan 2003, Ruszkowski & Begelman 2002, Fabian et al. 2001, Voigt et al. 2002, Gruzinov 
2002, Friaca 1986, Bertschinger & Meiksin 1986, Meiksin 1988, Bregman & David 1988). 

The role of thermal conduction in the ICM has been the subject of a long debate and, 
owing to the complex physics of MHD turbulence, the value of the effective conductivity 
remains uncertain. The thermal conductivity of an unmagnetised, fully ionised plasma was 
calculated by Spitzer (1962). Originally it has been thought that the magnetic field in 
clusters strongly supresses the thermal conductivity because the magnetic fields prevent an 
efficient transport perpendicular to the field lines. Even if the transport can be efficient along 
the magnetic field lines, the overall isotropic conductivity was thought to be many orders of 
magnitude less than the Spitzer value. This paradigm has been supported by a number of 
observations, such as sharp edges at so-called cold fronts, small-scale temperature variations 
in mergers and sharp boundaries around radio bubbles. It is thought that the existence of 
these sharp features precludes thermal conduction near the Spitzer value (Markevitch et al. 
2000, Vikhlinin, Markevitch & Murray 2001). 

Recent theoretical work by Narayan & Medvedev (2001), Malyshkin & Kulsrud (2001), 
Chandran et al. (1999), Chandran & Cowley (1998) and earlier work by Rechester & Rosen- 
bluth (1978) has shown that a turbulent magnetic field is not as efficient in suppressing 
thermal conduction as previously thought. It is argued that chaotic transverse motions of 
the tangled magnetic field lines can enhance the cross-field diffusion to an extent that the 
effective conductivity is of the order of the Spitzer value. Following this work, Zakamska & 
Narayan (2003) have looked for hydrostatically stable models in which the radiative cooling 
is exactly balanced by thermal conduction and where the temperatures on the inner and 
outer boundaries were fixed. For half of the clusters they investigated they found that a 
thermal conductivity of around 30 % of the Spitzer value yielded good fits to the observed 
profiles. However, for the other half of the clusters conduction alone appeared unlikely to 
halt a cooling catastrophe. Consequently, some form of heating is necessary such as me- 
chanical heating by AGN. The work by Zakamska & Narayan (2003) has motivated us to 
reexamine cluster models that takes into account thermal conduction as well as heating by 
central radio sources. We include a more generally valid treatment of radiative cooling by 
utilising a fit to the cooling function that takes into account line cooling. This becomes 
particularly important in the central regions of cooler clusters. 

Radio-loud AGN drive strong outflows in the form of jets that inflate bubbles or lobes. 
The lobes are filled with hot plasma, and can heat the cluster gas in various ways. The effect 
of hot bubbles on the ICM consists primarily of heating via pdV work and redistribution of 
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mass via buoyancy-driven mixing. Hydrodynamic simulations (Briiggen 2003) have shown 
that a significant fraction (~ 10 %) of the energy residing in radio lobes can be dissipated 
in the cluster gas. 

Clearly, the lifetime of the activity of the central AGN (as brief as ~ 10 4 — 10 5 yrs) is 
short compared to the evolutionary time scale of the cluster gas. Therefore, once the AGN 
stops supplying energy to the buoyant bubbles, the cluster gas will settle down once again 
and a full cooling flow may be re-established. The cooling gas flows to the centre of the cluster 
and may then trigger a further active phase of the AGN. Thus a self-regulating process with 
cooling periods alternating with brief bursts of AGN activity may be established (Quilis, 
Bower, & Balogh 2001, Voit & Bryan 2001). The rising bubbles uplift colder material from 
the vicinity of the AGN and thus disrupt the supply of fresh fuel for the radio jets. This 
feedback mechanism might automatically regulate the power of the radio jets. It has been 
found that 71 % of all cD galaxies at the centres of cooling-flow clusters show evidence of 
radio activity (Burns 1990). This fraction is higher than in non-cooling flow galaxies which 
again points towards the existence of some form of feedback mechanism where the cooling 
gas flows to the centre of the cluster and triggers a further active phase of the AGN. The 
bubbles rise at a speed that is comparable to the sound speed, which in turn is comparable 
to the dynamical speed. Therefore, the cooling timescale is much longer than the bubble 
rise timescale, and it is justifiable to treat feedback heating in a time- averaged sense. 

In this paper, we will present a simple model that incorporates, both, thermal conduction 
and heating by bubbles. We compute a class of models that are subject to radiative cooling, 
thermal conduction, heating by AGN and that are in hydrostatic equilibrium. We will 
investigate those parameters that yield plausible solutions and will compare our solutions to 
observations. Finally, we briefly discuss the origin of the energy that is conducted into the 
cluster. 



2. Model 

In our model of the cluster, we assume that the dark matter distribution is given by a 
modified NFW profile (Navarro et al. 1997) 

where r is the distance from the centre, r s the scale radius of the NFW profile and r c a soften- 
ing radius, within which the density becomes constant and which prevents the temperature 
from going to zero at r = 0. The core radius r c determines the shape of the potential near 
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the centre and here we adopt r c = r s /20 as recommended by Zakamska & Narayan (2003). 
We can express the characteristic mass M in terms of the commonly used concentration 
parameter c = (3M vil ./47r200p crit (-2)r^) 1 / 3 , where p crit is the critical density of the universe at 
the redshift of the cluster and M vir is the virial mass: 

n 2 / % /200\ c 3 
M = 2,r,W( Z ) (- j M1 + e) _ e/(1 + e) ■ (2) 

The conductive heat flux is given by 



F C = -V-(«VT), (3) 

where k is the thermal conductivity and T temperature. If thermal conduction is due to 
electrons, the conductivity according to Spitzer (1962) is given by 



k sp « 9.2 10~ 7 T 5/2 erg s" 1 K" 1 cm" 1 . (4) 

Here we seek an equilibrium model that is spherically symmetric, static, and time- independent. 
Moreover, we neglect magnetic fields and the self-gravity of the ICM. In spherical coordi- 
nates, the gravitational potential <3> is governed by the dark matter distribution 



1-f (r 2 $) = 4vrG PDM (r) = - , (5) 

r l dr v ' (r + r c )(r + r s ) z 

where G is the gravitational constant. Following Zakamska & Narayan (2003) the mass- 
temperature relation of Afshordi & Cen (2002) and the mass-scale relation of Maoz et al. 
(1997), can be used to determine Mo and r s for a given cluster. The equation of hydrostatic 
equilibrium can now be written as 



dp d$ 

where p is pressure and p gas density. Demanding that radiative cooling and heating are 
balanced by thermal conduction, we can write 



~^F c ) = - P C + n, (7) 

where pC and 7i are the volume cooling and heating functions, respectively, and F c is the 
conductive flux, which depends on the temperature gradient as stated above 
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«£ = , (8) 

where we write k as //%>, f being a suppression factor. We should point out that we omitted 
a term for the convective flux in equation (7). However, for the levels of heating considered 
here the ICM is convectively stable and the neglect of convection is a justifiable assumption. 

Pressure is related to n e and T via the ideal gas law 



p= pkT= l^ nek T, (9) 
pm u fx 



where m u is the atomic mass unit. 



For the volume cooling rate pC we use an approximation to the cooling function based 
on calculations by (Sutherland and Dopita 1993) 



pC = [C 1 (kT) a + C 2 {kTf + C 3 ]^n 2 e 10" 22 erg cm 3 s" 1 . (10) 

/x 

The units for kT are keV, n e is the electron number density, and fx and /x e denote the 
mean molecular weight per hydrogen atom and per electron, respectively. As in Zakamska 
& Narayan (2003) we use fx = 0.62 and /x e = 1.18, corresponding to a fully ionized gas 
with hydrogen fraction X = 0.7 and helium fraction Y = 0.28. For an average metallicity 
Z = 0.3Z & the constants are a = -1.7, /3 = 0.5, d = 8.6 x 10~ 3 , C 2 = 5.8 x 10~ 2 and 
C 3 = 6.4 x 10- 2 . 



The term 7i in equation (7) represents the energy input by a central AGN. The energy 
is deposited in the ICM by the rising bubbles and thus, averaged over time, the heating 
will be distributed in radius. To quantify this heating, we use a prescription proposed by 
Begelman (2001) which quantifies the heating of the ICM by the rising bubbles. Assuming 
that this heating mechanism reaches a quasi-steady state, the details of the bubble motions 
and geometry should cancel and the energy flux may be written as 



eocp b (r) (7b ^ 1)/7b , (11) 

where ph(r) is the partial pressure of buoyant gas inside bubbles at radius r and 7b is the 
adiabatic index of buoyant gas which we will take to be 4/3. Assuming that the partial 
pressure scales with the thermal pressure of the ICM, the volume heating function H can be 
expressed as 
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W~-Mr)V-— 5 = -MO f -77f. (12) 



47rr 2 \Po/ rdlnr' 

where p is the central pressure and h(r) is the normalization function 

h{ r) = J- 2 (l - e-^o )q -\ (13) 
where L is the luminosity of the AGN. The normalisation factor q is defined by 

P Y™- 1)/11 1 dlnp,, 



*=L \*) ^ (1 -^" ,)d "- (14) 

where r is the inner heating cutoff radius which is determined by the finite size of the cen- 
tral radio source. Here r was taken to be 10 kpc. This heating function mimics a possible 
feedback mechanism between the AGN and the ICM in the sense that the volume heating 
function does not depend on the physical conditions at the source alone but on the pressure 
gradient (Ruszkowski & Begelman 2002). Thus, it resembles thermal conduction, with the 
difference that the heating rate depends on the gradient of pressure rather than temperature. 



Equations (5) - (9) can be combined to yield a set of four first-order differential equations 
for n e , r 2 d§/dr, T and r 2 F c . We solve these equations as an initial- value problem subject to 
the initial conditions: n e (0) = n- in , r 2 d&/dr\ r=0 = 0, T(0) = T in and r 2 F c \ r=0 = 0, where n- in 
and T in are parameters. Other parameters in this simplified model are M , r s , f and L. We 
have investigated the dependence of the resulting density and temperature profiles on some 
of these parameters, which is discussed in the next session. 



3. Results and Discussion 



We integrate the system of four ordinary differential equations using a Runge-Kutta 
method and adopt a characteristic mass of Mq = 6.6 x 10 14 M Q and a scale radius of 
r s = 460 kpc (parameters inferred for the cluster Abell 1795, Zakamska & Narayan 2003). 
For a cluster with a central temperature of 2 keV and a central density of n e = 0.05 cm -3 the 
resulting density and temperature profiles are shown in Figures 1 and 2. Curves plotted in 
different styles correspond to different values of the suppression factor / and the luminosity of 
the central source L. The values assumed for L here correspond to typical energies supplied 
by the jet to the ICM of ~ 10 44 erg s" 1 (e.g. Owen, Eilek and Kassim 2000). In Fi gure 1 one 
can note that the temperature profile rises less steeply with radius if thermal conduction is 
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more efficient. In order to maintain a given central temperature, the temperatures gradient 
in the cluster has to be higher if conduction is more suppressed. The higher temperature 
gradient makes up for the smaller suppression factor so that enough energy is conducted 
inwards to balance the radiative losses. The highest value for the suppression factor that 
we have adopted here (0.5) is at the upper end of what is physically plausible. Even higher 
thermal conductivities seem very unlikely. Moreover, one can see, that the heating term has 
an effect similar to that of thermal conduction, in that a higher luminosity of the central 
source flattens out the temperature profile. Because the distributed heating depends on the 
pressure gradient, a higher luminosity of the central source can afford a smaller pressure 
gradient, leading to a shallower temperature distribution. Physically speaking, the smaller 
pressure gradient makes the energy transfer from the bubbles to the ICM less efficient. 

For higher values of / and L the electron number density is higher in the centre (< 100 
kpc) and drops off faster in the outer regions (see Fig. 2). 

Figures 3 and 4 are the corresponding figures for a cluster with a higher central temper- 
ature of 4 keV, all other parameters being the same. For this hotter cluster the luminosity of 
the central source has to be higher in order to produce the same relative effect as in the cooler 
cluster. In these plots we also included a model with a high luminosity of L = 3 x 10 45 erg 
s _1 . In this case the temperature distribution is nearly isothermal. 

Clearly, by including a physically motivated heating term, there are more degrees of 
freedom to fit the density and temperature profiles inferred from observations. This can be 
demonstrated at the example of Hydra A which is a well-studied cooling flow cluster with a 
FR I radio source (3C218) at its centre. Zakamska & Narayan (2003) have shown that Hydra 
A cannot be fitted with a conduction model alone unless one allows for unplausibly high 
values for k. The NFW parameters for Hydra A are r s = 370 kpc and M = 2.9 x 10 14 M & . 
The central electron density is n e (0) = 0.08 cm' 3 , the temperature at the centre is T in = 
3.11 keV and at 190 kpc is T(190 kpc) = 4.04 keV. In Figure 5 we show the temperature 
profile as inferred from CHANDRA observations (David et al. 2001) together with our best 
fit for / = 0.5 (L = 4.0 x 10 45 erg s~ 4 ). As is readily seen, the inclusion of the heating term 
produces a good fit to the data without having to invoke an abnormally high conductivity. 

A question that we have not addressed so far is the question of the origin of the energy 
that is conducted inwards in order to keep the cluster from collapsing. If we take a cluster 
model with no heating, i.e. L = 0, and a suppression factor of / = 0.1 (all other parameters 
being those of Fig. 1), the asymptotic value of r 2 F c at large radii is 6.2 x 10 44 erg s _1 . This 
means that over a period of 5 Gyrs a total energy of ~ 10 62 erg has to be conducted into the 
cluster which amounts to ~ 5 % of the total thermal energy of the cluster (3MokT /2fim p ). 
For higher values of / this percentage increases. Ultimately, this energy must be provided 
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by infalling gas at the accretion shock. 

It was pointed out by Loeb (2002) that a significant thermal conductivity will also 
lead to energy transport out of the cluster because beyond the temperature maximum of 
the cluster the temperature begins to drop again and the temperature gradient is reversed. 
Therefore, conduction is also reversed and heat is transported from the cluster to its sur- 
rounding envelope in these regions. As estimated by Loeb (2002), heat conduction must be 
suppressed significantly over the Spitzer value in the outer cluster regions or else the cores 
of X-ray clusters would have cooled significantly over their life times. 

3.1. Stability of the solutions 

We have not performed a stability analysis of the model presented here. Instead we 
refer the reader to relevant works that have been published by others. The thermal stability 
of stratified atmospheres in the presence of radiative cooling and heat conduction was first 
analysed by Field (1965) and later generalised by Balbus (1988) and Balbus and Soker (1989) 
who performed a Lagrangian stability analysis. Recently, the issue of thermal instability in 
clusters was revisited by Kim & Narayan (2003) who examined the global stability of the 
models found by Zakamska & Narayan (2003). On the basis of a Lagrangian perturbation 
analysis they concluded that the growth time of the most unstable global radial mode was 
~ 6 - 9 times longer than the growth time of the local isobaric modes at the centre if 
the conductivity was 20 % to 40 % of the Spitzer value. Typical growth times are thus 
comparable with the time since the last cluster merger. Hence it is argued that the cluster 
is thermally stable if there is a sufficient amount of thermal conduction. The presence of a 
heating term as implemented in our model will act to increase the stability of the cluster. 

3.2. Summary 

We have devised a one-dimensional hydrostatic model for the ICM where radiative losses 
are balanced by, both, thermal conduction and heating by a central source. By including 
a physically motivated heating term we have shown that it is possible to fit cluster profiles 
without having to invoke unplausibly high values for the thermal conductivity. 

We should concede, however, that our model is very simplified: It assumes that the 
cluster is in hydrostatic equilibrium and spherically symmetric, it does not allow the gas 
to condense and drop out of the flow, it ignores magnetic fields and includes a heuristic 
treatment of heating by AGN. But despite its simplicity, only a few free parameters suffice 
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to reproduce the observed profiles of X-ray clusters. 



I thank Simon White, Eugene Churazov, Bill Forman and Christian Kaiser for helpful 
discussions. 
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Fig. 1. — Temperature profiles for a cluster with central temperature T in = 2 keV and central 
density of n e = 0.05 cm" 3 with / = 0.1, L = 3 x 10 43 erg s" 1 (solid) and / = 0.1, L = 3 x 10 44 
erg s" 1 (dotted), / = 0.5, L = 3 x 10 43 erg s" 1 (dashed) and / = 0.5, L = 3 x 10 44 erg s" 1 
(dash- dot). 
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Fig. 2. — Electron number density profiles corresponding to the models presented in Fig. 1. 
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Fig. 3. — Temperature profiles for a cluster with central temperature Ti n = 4 keV and central 
density of n e = 0.05 cm" 3 with / = 0.1, L = 3x 10 43 erg s" 1 (solid) and / = 0.1, L = 3 x 10 44 
erg s" 1 (dotted), / = 0.1, L = 3 x 10 45 erg s" 1 (dashed), / = 0.5, L = 3 x 10 43 erg s" 1 
(dash-dot) and / = 0.5, L = 3 x 10 44 erg s -1 (dash- dot- dot). 
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Fig. 4. — Electron number density profiles corresponding to the models presented in Fig. 3. 
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Fig. 5. — Gas temperature profile of the Hydra A cluster. The crosses show the results of 
ACIS data published in David et al. (2001). The line corresponds to the fit with / = 0.5 
and L = 4.0 x 10 45 erg s _1 . 
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